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Abstract. A spherically symmetric model is presented 
for the interaction of a pulsar wind with the associated 
supernova remnant. This results in a pulsar wind nebula 
whose evolution is coupled to the evolution of the sur- 
rounding supernova remnant. This evolution can be di- 
vided in three stages. The first stage is characterised by 
a supersonic expansion of the pulsar wind nebula into the 
freely expanding ejecta of the progenitor star. In the next 
stage the pulsar wind nebula is not steady; the pulsar wind 
nebula oscillates between contraction and expansion due 
to interaction with the reverse shock of the supernova rem- 
nant: reverberations which propagate forward and back- 
ward in the remnant. After the reverberations of the re- 
verse shock have almost completely vanished and the su- 
pernova remnant has relaxed to a Sedov solution, the ex- 
pansion of the pulsar wind nebula proceeds subsonically. 
In this paper we present results from hydrodynamical sim- 
ulations of a pulsar wind nebula through all these stages 
in its evolution. The simulations were carried out with the 
Versatile Advection Code. 



1. Introduction 

The explosion of a massive star at the end of its life as 
a supernova releases an amount of energy roughly equal 
to 10^'^ erg. Roughly 99% of this energy is radiated away 
in the form of neutrinos as a result of the deleptonization 
of the ~ lAf© stellar core as it collapses into a neutron 
star. The remaining (mechanical) energy of about lO'^^ erg 
is contained in the strong shock propagating through the 
stellar mantle, and ultimately drives the expansion of the 
supernova remnant (SNR). 

In those cases where a rapidly rotating neutron star 
(pulsar) remains as a 'fossil' of the exploded star, a pulsar 
wind, driven by the spindown luminosity of the pulsar, can 
be formed. The precise magnetospheric physics leading to 
such a pulsar wind is not fully understood, but it is be- 
lieved that a major fraction of the spin-down luminosity 
of the pulsar is converted into the mechanical luminosity 
of such a wind. 



The total rotational energy released by a Crab- like pul- 
sar over its lifetime is of order 10^^ — 10^° erg. This is much 
less than the mechanical energy of ~ 10^^ erg driving the 
expansion of the SNR. Therefore, the dynamical influence 
of the pulsar wind on the global evolution of the supernova 
remnant itself will be small. From an observational point 
of view, however, the presence of a pulsar wind can lead to 
a plerionic supernova remnant, where the emission at ra- 
dio wavelengths shows an extended, flat-spectrum central 
source associated with a Pulsar Wind Nebula (PWN) . The 
best-known example of such a system is the Crab Nebula, 
and about a half-dozen other PWNs are known unam- 
biguously around pulsars from radio surveys (e.g. Frail & 
Scharringhausen, 1997). These surveys suggest that only 
young pulsars with a high spindown luminosity produce 
observable PWNs at radio wavelengths. At other than ra- 
dio wavelengths, in particular X-rays, there are about ten 
detections of PWN around pulsars both in our own galaxy 
and in the large Magellanic cloud (LMC) (e.g. Helfand, 
1998; Table 1 of Chevaher, 2000). 

The expansion of an isolated SNR into the general 
interstellar medium (ISM) can be divided in four differ- 
ent stages (Woltjer, 1972; see also CiofS, 1990 for a re- 
view): the free expansion stage, the Sedov- Taylor stage, 
the pressure-driven snowplow stage and the momentum- 
conserving stage. In the models presented here wc will only 
focus on the evolution of a pulsar wind nebula (PWN) 
during the first two stages of the SNR: the free expan- 
sion stage and the Sedov- Taylor stage. We will assume 
that the pulsar is stationary at the center of the remnant, 
excluding such cases as CTB80 (e.g. Strom, 1987; Hester 
& Kulkarni, 1988), PSR1643-43 and PSR1706-44 (Frail et 
al., 1994), where the pulsar position is significantly excen- 
tric with respect to the SNR, presumably due to a large 
kick velocity of the pulsar incurred at its birth, assuming 
of course that SNR and pulsar are associated and we are 
not dealing with a chance alignment of unrelated sources. 
The case of a pulsar moving through the remnant with a 
significant velocity will be treated in a later paper. 

In this paper we compare (approximate) analytical ex- 
pressions for the expansion of a PWN in a supernova rem- 
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nant with hydrodynamical simulations carried out with 
the Versatile Advection Code |^ (VAC). We confirm earlier 
analytical results (Reynolds & Chevalier, 1984; Chevaher 
&; Fransson, 1992) which state that the PWN is expand- 
ing supersonically when it is moving through the freely 
expanding ejecta of the SNR. Due to deceleration of the 
expanding SNR ejecta by the interstellar medium (ISM), 
a reverse shock propagates back to the center of the SNR 
(e.g. McKee, 1974 Cioffi et al. 1988)). Due to the presence 
of reverberations of the reverse shock in the SNR , the 
expansion of the PWN goes through an unsteady phase 
when this reverse shock hits the edge of the PNW. Af- 
ter these reverberations have decayed, the expansion of 
the PWN through the ejecta of the SNR progenitor star 
continues subsonically with the PWN almost in pressure 
equilibrium with the interior of the SNR. 

This paper is organised as follows. In sections 2 
and 3 we discuss the aforementioned two stages of the 
PWN/SNR system. In section 4 the hydrodynamical sim- 
ulations will be presented and compared with the analyt- 
ical expressions from section 2 and 3. 

2. Pulsar Wind Nebula in a freely expanding 
Supernova Remnant 

In the early stage of the evolution of a PWN, the SNR 
consists mostly of the stellar ejecta expanding freely into 
the interstellar medium. The PWN expands into these 
ejecta. The sound velocity in the interior of the SNR is 
much smaller than the expansion velocity of the PWN. 
The supersonic expansion of the PWN results in a shock 
propagating into the ejecta (see figure 1). 

An analytical equation for the radius of this shock can 
be derived for a constant spindown luminosity. Using this 
solution, the assumption of supersonic expansion will be 
checked a posteriori. For simplicity we assume that the 
ejecta have a uniform density, 
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We assume that the stellar ejecta swept up by the strong 
shock which bounds the PWN collect in a thin shell, and 
that this material moves with the post-shock velocity. Ne- 
glecting the contribution of the thermal energy we can 
write the total (kinetic) energy of this shell, iSsheii, as: 
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where 



Msw(i) = Mo, 



Rpwn 



(6) 



is the ejecta mass swept up by the pulsar wind nebula. 
In deriving the post-shock velocity, we assumed that the 
ejecta behave as an ideal non-relativistic gas with adia- 
batic heat ratio 7ej =5/3 and used the Rankine-Hugoniot 
jump conditions for a strong shock. 

The interior of the PWN is dominated by thermal en- 
ergy. The sound speed in a realistic PWN is close to the 
speed of light c, while the expansion velocity is much less 
than c. Perturbations in the pressure will be smoothed out 



^ -Rpwn/c, much less 
-^pwn /-^pwn • There- 
fore, we can assume a nearly uniform pressure Ppwn in the 
PWN. The internal energy of the PWN then equals 
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than the expansion time scale texp ' 
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and a linear velocity profile as a function of radius. 



Here we take 7pwn = 4/3 because the pulsar wind nebula 
material is relativistically hot. The pressure of the inte- 
rior of the PWN must roughly equal the pressure in the 
shocked ejecta just downstream of the outer shock of the 
pulsar wind nebula at Rpwn- 
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with i?oj = Vot the radius of the front of the ejecta. The 
value of Vq is determined by the requirement that the 
kinetic energy of the ejecta equal the total mechanical 
energy Eq of the SNR: 
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Combining these relations yields: 
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^ See http://www.phys.uu.nl/ toth/ 



(^3) Energy conservation for the PWN system reads: 

S,hcll(<) + Epwnit) = Einitit) + Lot . (10) 
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Here i?init(i) is the kinetic energy which the swept-up 
ejecta would have if they were freely expanding. This 
quantity can be obtained by integrating the kinetic en- 
ergy density of ejecta in a sphere with radius r < i?pwn if 
there was no PWN. This yields: 

i?i„itW = i?o (^)'- (11) 

After some algebra using the equations (|l|)-(pT|) one can 
obtain a power-law solution for the radius of the pulsar 
wind bubble: 

i?pw„(t) - C (^^^ ^ Vot cx t^/s, (12) 
where C is a numerical constant of order unity: 

5 4 

with 70j=g, 7pWn = g 

Reynolds and Chevalier(1984) already obtained this 
expansion law. It can easily be checked that 
the expansion velocity obtained in this manner is indeed 
much larger than the sound velocity in the freely expand- 
ing supernova remnant. 

3. Pulsar Wind Nebula in a Sedov- Taylor 
remnant 

3.1. Pulsar Wind Nebula expansion for a constant wind 
luminosity 

Towards the end of the free expansion stage a reverse 
shock is driven deep into the interior of the SNR. This 
reverse shock reheats the stellar ejecta, and as a result 
the sound velocity increases by a large factor. When the 
reverberations due to reflections of the reverse shock have 
almost completely dissipated, one can approximate the in- 
t erior of the SNR by using the analytical Sedov solution 
( ^edov, 1958]) . 

The interaction with the reverse shock influences the 
evolution of the pulsar wind nebula quite dramatically. 
CiofH et al. (1988) have already shown in their ID simula- 
tion of a pure shell SNR that the reverse shock gives rise to 
all kinds of sound waves and weak shocks traveling back 
and forth through the ejecta before the interior relaxes 
towards a Sedov solution. We will show that during the 
process of relaxation the radius of the pulsar wind nebula 
contracts and expands due to reverberations of the reverse 
shock. Compression waves are partly reflected and partly 
transmitted at the edge of the PWN. We will come back 




Fig. 1. Schematic representation of PWN in a freely ex- 
panding SNR. There are a total of four shocks and two 
contact discontinuities. From left to right one can see: 
the pulsar wind termination shock i?ts (dashed line), the 
first contact discontinuity i?cdi (dotted line) separating 
shocked pulsar wind material from shocked ejecta, the 
PWN shock i?pwn (solid line) bounding the PWN. For 
the SNR we have a reverse shock i?,,, (dashed line), the 
second contact discontinuity i?cd2 (dotted line) separating 
shocked ejecta from shocked ISM, and the SNR shock Rgnr 
(solid line) which is the outer boundary of the PWN/SNR 
system. 



to this point when we discuss results from hydrodynam- 
ics simulations in section 4, which allow a more detailed 
picture of this process. 

In this Section we consider a fully relaxed Sedov SNR. 
The PWN expands subsonically into the remnant because 
the interior of the SNR has been re-heated by the reverse 
shock. For the case of a constant (mechanical) luminosity 
driving the pulsar wind an analytical expression for the 
radius of the PWN can be easily obtained. In this stage 
of the PWN evolution, we associate its radius -Rpwn with 
the contact discontinuity separating pulsar wind material 
from the ejecta of the progenitor star (see figure 2). 

We first present an order-of-magnitude calculation 
which leads to the correct power-law solution for the ra- 
dius of the PWN. The assumption of subsonic expan- 
sion implies approximate pressure equilibrium between the 
wind material and the stellar ejecta at the edge the PWN. 
In the interior of the SNR the pressure scales as 

Psnr « Eo/R^ . (14) 
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Here Eth is the thermal energy of the PWN, Pi its internal 
pressure, and Vpwn its volume. This yields the following 
equation describing the energy balance of a slowly expand- 
ing PWN: 



d / 47r Pi-Rpwn 



dt \ 3 (7, 
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or equivalently 



47r 7pwn-Pi-Rp 
3 (7pwn 
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This equation has a power-law solution for i?pwn(i) pro- 
vided the internal pressure Pi{t) in the SNR behaves as a 
power-law in time so that the relation 



Fig. 2. Schematic representation of a PWN in a Sedov 
SNR. There are a total of 2 shocks and 2 contact dis- 
continuities. From left to right one can see: the pulsar 
wind termination shock i?ts (dashed line), the first con- 
tact discontinuity i?pwn (dotted line) separating shocked 
pulsar wind material from shocked ejecta, bounding the 
PWN. Furthermore there is another contact discontinuity 
i?cd (dotted line) separating shocked ejecta from shocked 
ISM, and the SNR shock i?sni. (solid line) which bounds 
the PWN/SNR system. 



On the other hand, the pressure in the interior of the PWN 
scales as 



Ppwn oc Lot/Rf 
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can be satisfied. For a Sedov SNR expanding into a uni- 
form ISM one has P\ cx and one finds: 
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If i?pwn ^ ^snr we Can use the central pressure from the 
Sedov solution with yism = 5/3 for the interior pressure in 
the SNR which confines the PWN (e.g. Shu, 1992): 



with Lq the mechanical luminosity driving the wind. Pres- 
sure equilibrium at the contact discontinuity at i?pwn im- 
plies the following relation for the radius of the PWN as 
a function of time: 

i?pwn(0 = C (^^^ ^ RsnAt) (X (16) 

with the constant of proportionality C to be determined 
below. 

A more detailed derivation uses the first law of ther- 
modynamics, assuming once again a constant energy input 
Lq into the PWN by the pulsar-driven wind: 

dEth = Lodt- Pi dVpwn- (17) 



m) = Ps„r(i) ^ 0.074 (^-^^ oc . (23) 

We find the same result for Ppwn jt) as in the order-of- 
magnitude calculation (Eqn. Hq), determining the con- 
stant in that expression as C* ~ 0.954 for a non-relativistic 
fluid (7pwn = 5/3) and C ~ 0.851 for a relativistic fluid 
(7pwn = 4/3). By comparing the sound speed with the ex- 
pansion velocity at the edge of the PWN, we confirm that 
the expansion remains subsonic. 

An alternative derivation of the PWN expansion law 
uses the Kennel-Coroniti model for a highly relativistic 
pulsar-driven pair wind. This wind is terminated by a 
strong MHD shock which decelerates the fluid to a non- 
relativistic expansion speed (Rees & Gunn, 1974). Kennel 
& Coroniti (1984, hereafter K&C) constructed a steady. 
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spherically symmetric MHD model for the Crab nebula 
which includes these characteristics. We use their model 
in the hydrodynamical limit by considering the case 



Poynting flux 







particle energy flux 
K&C assume a constant wind luminosity, 
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where ni is the proper density just in front of the termi- 
nation shock, u — Ti V ~ Tic is the radial four-speed of 
the wind and i?ts is the distance from the pulsar to the 
termination shock. Because the wind is assumed to con- 
sist solely of a positronic plasma, m is the electron mass. 
The pulsar wind is highly relativistic (Fi :s> 1) and the 
thermal and rest energy of the particles can be neglected 
compared with the bulk kinetic energy. The total number 
of particles emitted into the PWN then equals: 



N{t) 



Lot 



(26) 



It is believed that the bulk Lorentz factor Fi « 10^, but 
we will see that for purposes of the PWN evolution its 
precise value is not important, because it cancels in the 
final result. 

After the termination of the cold wind by a strong stand- 
ing shock at some radius i?ts, the wind flow is subsonic, 
with a sound speed close to the speed of light. Assum- 
ing that the shock radius i?ts is much smaller than the 
radius -Rpwn of the PWN, and assuming a uniform den- 
sity 71.2 and uniform pressure P2 inside the PWN, particle 
conservation implies 



Fig. 3. Comparison between results from numerical simu- 
lations and analytical result for the radius of the PWN, i.e. 
equation ( p^ . The dashed line indicates the radius for the 
PWN obtained from numerical simulations. The solid line 
corresponds to equation ( |T2| ) with C ~ 0.941, as appropri- 
ate for 7pwn = 5/3. The different physical parameters are 
as indicated in Table 1 (Simulation 1). The injected mass 
of the pulsar wind has been chosen in such a way that the 
termination velocity of the pulsar wind equals the speed of 
light. One can see that in the simulation the radius i?pwn 
is about 10 % smaller than predicted by the analytical re- 
sult, but the power-law behaviour i?pwn oc t^^^ is correctly 
reproduced. 



the downstream pressure P2. At the outer edge of the 
PWN this pressure must approximately equal the pres- 
sure Psnr(i) at the center of the SNR as given by Eqn. 
(p3|). Using ( p8| ) and ( p9| ) the inner and outer boundary 
conditions imply the following relation for the termination 
shock radius: 



— n2 Rl^a[t) - 



Fi mc^ 



(27) 



From the K&C model we take the following relationships, 
valid at the strong relativistic termination shock at the 
inner edge of the PWN in the hydrodynamical limit: 



2 2 
mc 
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The subscripts 1 and 2 label upstream and downstream 
parameters on either side of the termination shock. Using 
these jump conditions together with equations (|2^) and 
(p6|) we can express Rpwnit) as a function of i?ts(i)- 
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The pressure inside the PWN is nearly uniform. At the 
termination shock (inner edge of the PWN) it must equal 



Rts{t) ~ 0.847 



Lq 
Eoc 



1/2 



(30) 



It is now straightforward to obtain the radius of the PWN 
from (^ and ( |30|) . The resulting expression for i?pwn sat- 
isfies equation ( |16[ ) with C ~ 0.911. This derivation based 
on (ram) pressure balance at the inner and outer edges of 
the pulsar wind nebula confirms our earlier result obtained 
from overall energy conservation. 



3.2. Pulsar Wind Nebula expansion for varying wind 
luminosity 

The constant wind luminosity assumption is not very re- 
alistic by the time the effects of the reverse shock and its 
associated reverberations have vanished. The spin-down 
luminosity of the pulsar is more realistically described by 
the luminosity evolution from a rotating magnetic dipole 
model: 
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L{t) = ^\ 2 - (31) 



Therefore we now consider the more realistic case of a 
time-dependent luminosity given by The energy bal- 
ance equation for the PWN reads: 

_d_ / 47r PjRl^n \ _ Lq _ 

dtys (7pwn - i)y C^ + t/rf 

AnRl^^ P, . (32) 

We solve this equation numerically using a fourth-order 
Runge-Kutta method (e.g. Press et al., 1992). As an initial 
condition we take the radius of the PWN equal to zero at 
the start of the evolution, neglecting the initial stage when 
the PWN is expanding supersonically. For the pressure Pi, 
we use the pressure at the center of the Sedov SNR (p3|). 
We find that the solution for i?pwn converges to i?pwn oc 
i"'^ on a time scale much larger than the typical time scale 
for the reverse shock to hit the edge of the PWN. Figure 
9 shows this semi-analytical result together with results 
from hydrodynamical simulations. For the semi- analytical 
equation we use 7pwn = 5/3, because the hydrodynamics 
code also uses this value (see section 4.1 below). 

4. Numerical simulations 
4.1. Method 

Our simulations were performed using the Versatile Ad- 
vection Code (VAC, Toth 1996) which can integrate the 
equations of gas dynamics in a conservative form in 1, 2 or 
3 dimensions. We used the TVD-MUSCL scheme with a 
Roe-type approximate Riemann solver from the numerical 
algorithms available in VAC (Toth and Odstrcil, 1996); a 
discussion of this and other schemes for numerical hydro- 
dynamics can be found in LeVeque (1998). In this paper 
our calculations are limited to spherically symmetric flows. 
We use a uniform grid with a grid spacing chosen sufli- 
ciently fine to resolve both the shocks inside the PWN 
and the larger-scale shocks associated with the SNR. Ta- 
ble 1 gives the physical scale associated with the grid size 
for the simulations presented here. An expanding SNR is 
created by impulsively releasing the mechanical energy of 
the SN explosion in the first few grid cells. The thermal 
energy and mass deposited there lead to freely expanding 
ejecta with a nearly uniform density, and a linear velocity 
profile as a function of radius. 

A realistic shocked pulsar wind is presumably highly rel- 
ativistic with an adiabatic heat ratio 7pwn = 4/3. The 
(shocked) stellar ejecta on the other hand are non - rela- 
tivistic with 7oj = 5/3. 

The VAC code does not currently include relativistic 
hydrodynamics. Therefore, the best approach available to 
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Fig. 4. Pressure profile of the PWN/SNR system as a 
function of radius, at time t = 1000 years after the SN 
explosion. Physical parameters are as indicated in Table 1 
(Simulation 2). Moving outwards in radius one can see the 
wind termination shock, the shock bounding the PWN, 
the reverse shock of the SNR and the shock bounding the 
SNR. The interior of the PWN is nearly isobaric. There 
is a sudden increase in pressure of the ejecta behind the 
SNR reverse shock. 
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Fig. 5. Density profile for the same PWN/SNR system as 
in figure 4. 



us is to keep 7pwn = 5/3, but to take a luminosity for 
the pulsar wind, L(t), and an associated mass injection, 
Mpvj(t), such that the terminal velocity obtained from 
these two parameters, 

voo = y^2L(i)/Mcj(0 , (33) 

roughly equals the speed of light. Since the pulsar wind 
material downstream of the termination shock moves with 
only a mildly relativistic bulk speed we expect our re- 
sults to be qualitatively correct. Thermal energy and mass 
are deposited continuously in a small volume as a source 
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Fig. 6. Sound velocity profile as a function of radius, 
for the same case as in figures 4 and 5. Because of the 
Rankine-Hugoniot jump conditions at the wind termina- 
tion shock, the sound velocity of the shocked wind ma- 
terial in the PWN bubble is close to the speed of light. 
Behind the contact discontinuity, where the bubble con- 
sists of swept-up ejecta, the sound speed has a smaller 
value. 



for the wind. The hydrodynamics code then develops a 
steady wind reaching the terminal velocity Voo well before 
the (cold) wind is terminated by the standing termination 
shock. 

We trace the total mass injected into the PWN by the 
pulsar wind in order to determine the radius of the con- 
tact discontinuity which separates the pulsar wind mate- 
rial from the SN ejecta (i?cdi in figure 1 and i?pwn in figure 
2). We also determine the position of the shock bounding 
the PWN during the stage of supersonic expansion. This 
enables us to compare the numerical results with the ana- 
lytical expressions derived in sections 2 and 3 for the PWN 
radius. 

As a test of the code we have calculated a pulsar wind 
driven by a constant luminosity Lq (Simulation 1 in Table 
1). We let the PWN evolve until the reverse shock propa- 
gating in the SNR hits its outer edge. Figure 3 shows the 
radius of the shock of the PWN together with the analyt- 
ical equation (p^). We take 7pwn = 5/3 in the analytical 
expressions for comparison with the numerical results. Al- 
though the analytical result of Eqn. (^6|) is not reproduced 
exactly (the radius is about 10% smaller), the power-law 
expansion law i?pwn t^/^ is reproduced. As we will show 
in section 4.2, the pressure inside the bubble is larger than 
the one used to derive the equation, explaining the differ- 
ence between the analytical and numerical results. 
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Fig. 7. Velocity profile for the PWN/SNR system with the 
same parameters. The terminal wind velocity, Voc, is close 
to the speed of light. The large jump in the velocity at a 
radius ~ 4 pc is the reverse shock which is still propagating 
forwards in the laboratory frame. The velocity jump of the 
PWN shock at radius ^ 2.5 pc is much smaller. 



4-.2. Evolution of the PWN-SNR system into the Sedov 
phase 

Our simulations of the evolution of a pulsar wind nebula 
inside a supernova remnant employ the parameters listed 
in Table 1. 

In the early stage of its evolution the PWN is bounded 
by a strong shock propagating through the ejecta of the 
progenitor star. In figures 4-7 one can clearly identify the 
four shocks indicated schematically in figure 1. Moving 
outward in radius one first encounters the pulsar wind ter- 
mination shock; this termination shock is followed by the 
PWN shock. In the sound velocity profile of figure 6 one 
can see a large jump between these two shocks: the con- 
tact discontinuity separating shocked pulsar wind material 
from shocked ejecta. Further outward one encounters the 
SNR reverse shock, which at this stage of the SNR evolu- 
tion is still moving outwards from the point of view of a 
stationary outside observer. The whole PWN-SNR system 
is bounded by the SNR blast wave. 

Figure 9 shows the evolution of the contact discontinu- 
ity radius i?cd, which can be identified with the radius of 
the PWN in the subsonic expansion stage. One can clearly 
see the moment at t ~ 1.75 kyr when the reverse shock 
hits the edge of the PWN: the expansion becomes un- 
steady with the PWN contracting and expanding due to 
the interaction with the pressure pulses associated with 
the reverberations of the reverse shock. When these re- 
verberations have almost dissipated the expansion of the 
PWN relaxes to a steady subsonic expansion. In this stage, 
we can fit the radius of the PWN obtained from the simu- 
lations with the (semi-) analytical solution obtained from 
a numerical integration of equation ( |l8|) , as shown in this 
figure. 
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Table 1: Simulation parameters 





Simulation 1 


Simulation 2 


Simulation 3 


Explosion energy £"0 (erg) 


10-51 


10^1 


10-54 


Ejecta mass Mej(M0) 


3 


3 


3 


Pulsar wind luminosity Lq (erg/s) 


5 X 10^*^ 


1038 


5 X 10^*^ 


Spin-down time r (yr) 


00 


600 


600 


ISM mass density po (g/cm"^) 


10-24 


10-24 


10-24 


Number of grid cells 


5000 


5000 


3000 


Grid size (pc) 


0.002 


0.002 


0.01 




Fig. 8. The radius of the PWN contact discontinuity as a 
function of time (sohd line) . We compare with the semi- 
analytical solution from equation (|3^) (dashed line) . Here 
one can see that the expansion of the PWN is unsteady, 
due to the reverberations of the reverse shock. This sim- 
ulation was done with the parameters listed in Table 1 
(Simulation 3). 



The interaction of the PWN with the reverse shock and 
the associated reverberations is quite complicated. We will 
therefore describe this process in more detail. 

4-.3. The influence of reverse-shock reverberations 

The reverse shock initially encounters the PWN in its 
supersonic expansion stage. After the collision between 
the reverse shock and the PWN shock a reflected shock 
propagates back towards the outer (Sedov- Taylor) blast 
wave of the SNR. A transmitted shock propagates into 




4 6 
Time (in 1 0^ years) 
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Fig. 9. The radius of the PWN as a function of time, 
together with the ratio of the total energy in the PWN 
bubble with respect to the total energy input by the pulsar 
wind. The solid line represents the radius of the contact 
discontinuity of the PWN (in arbitrary units), the open 
squares represent the aforementioned ratio of energy. This 
simulation was done with the parameters as listed in Table 
1 (Simulation 3). 



the shocked ejecta inside the PWN. When this shock hits 
the contact discontinuity bounding the pulsar wind mate- 
rial a similar reflection/transmission event occurs: a shock 
moves radially outwards, and a compression wave moves 
into the pulsar wind material. The latter wave is rapidly 
dissipated in the pulsar wind bubble because of the high 
sound speed in the shocked pulsar wind. After a few sound 
crossing times the pulsar wind bubble contracts adiabati- 
cally in response to the pressure increase inside the SNR. 
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After this contraction it regains pressure equilibrium with 
the surrounding SNR and the PWN expands subsonically 
henceforth. This chain of events can be clearly seen in fig- 
ure 8 where we plot the radius of the PWN. The whole 
process takes a time comparable with the duration of the 
initial supersonic expansion stage. 

4-.4- Subsonic expansion stage 

When the PWN has more or less relaxed to a steady sub- 
sonic expansion the PWN has gained energy as a result 
of the interaction with the reverse shock. Consequently, 
the radius of the PWN is roughly 20% larger than the 
value predicted by the semi-analytical solution obtained 
from Eqn. (|l8|) in Section 3.2. In figure 9 we show the 
ratio between the (mostly thermal) energy of the pulsar 
wind bubble, i.e. the part of the PWN that consists of 
shocked pulsar wind material, and the total mechanical 
energy deposited by the pulsar. One can clearly see the 
increase in the energy content of the pulsar wind bubble. 
A large fraction of the energy deposited by the pulsar wind 
in the stage when the expansion is supersonic is contained 
in the kinetic energy of the shocked stellar ejecta in the 
PWN shell. When the reverse SNR shock is interacting 
with the PWN bubble, energy is apparently transferred 
from this thin shell to the interior of the bubble through 
the dissipation of the waves transmitted into the bubble. 

5. Conclusions and discussion 

We have considered a spherically symmetric PWN/SNR 
system in the early and middle stages of its evolution, 
well before cooling of the SNR shell becomes dynamically 
important and before a significant disruption of spherical 
symmetry due to a possible (large) kick velocity of the 
pulsar can take place. The expansion of the PWN is cou- 
pled with the dynamics of the expanding SNR, leading to 
two distinct evolutionary stages separated by an unsteady 
transition phase: 

— When the PWN is surrounded by the freely expanding 
ejecta of the SNR, the expansion of the PWN is super- 
sonic. In this stage the pressure in the interior of the 
PWN bubble is slightly larger than one would expect 
from ram pressure of the surrounding ejecta alone, us- 
ing the Rankine-Hugoniot relations at the PWN shock. 
This is due to the thin shell of shocked, swept-up ejecta 
which needs to be accelerated by the outward force due 
to the interior pressure of the PWN. 

— This stage of supersonic expansion is ultimately fol- 
lowed by a subsonic expansion of the PWN. This hap- 
pens after the reverse shock has encountered the shock 
of the PWN. 

— The transition between these two stages is unsteady 
due to the interaction of the PWN with the reverse 
shock and its associated reverberations. From the hy- 
drodynamical simulations we see that the time scale for 



adjustment to the pressure of the surrounding SNR is 
determined by the sound speed of the ejecta shocked 
by the PWN in the first stage. 

Two of the prototypes of plerionic SNRs are the Crab 
Nebula and 3C58. In the Crab, there is a decrease in the 
radio flux of 0.167 ± 0.015 % yr'^ (Aller & Reynolds 
1985). By contrast, 3C58 shows an increase in its flux den- 
sity at radio frequencies between 1967 and 1986 (Green 
1987). This increase might be the result of the reverse 
shock which has encountered the PWN shock around 
3C58; the PWN is being compressed and therefore the 
flux density is going up. 

Our numerical simulations are different from the re- 
sults presented by Jun (1998). This author concentrates on 
the details of the PWN in the supersonic expansion stage, 
and in particular on the formation of Rayleigh-Taylor fin- 
gers in his two-dimensional simulations. Our simulations 
include the whole supernova remnant, but can not ad- 
dress the development of Rayleigh-Taylor instabilities due 
to our assumption of spherical symmetry. 

In future work we will discuss how these results change 
when the influence of a significant kick velocity of the pul- 
sar is taken into account. If this is taken into account, 
the model presented here will lose its validity at a cer- 
tain time: one can calculate when the motion of a pulsar 
will become supersonic in a Sedov stage. One can show 
that this will happen when the pulsar is about halfway 
from the explosion center to the edge of the SNR: a bow 
shock is expected to result from this and clearly the model 
presented here will break down. Observationally there is 
evidence that this is the case for the pulsar associated with 
the SNR W44. 
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